##################################################
## Author: Stephanie L. DeMora (1), Jennifer L. Merolla (1), Brian Newman (2), and Elizabeth J. Zechmeister (3)
## Project: "Reducing Mask Resistance with Value Consistent Messages: A Study of White Evangelical Christians"
## Purpose: Replication File for Figures and Tables
## Date: Feb 19th 2021
##################################################


# Read in packages:
options(scipen=999)
library(haven)
library(ggplot2)
library(broom)
library(knitr)
library(kableExtra)
library(dotwhisker)
library(sjPlot)
library(sjlabelled)
library(Rmisc)
library(stringr)
library(tidyverse)
library(questionr)
library(weights)
library(jtools)
library(margins)
library(interplot)
library(nnet)
library(ggpubr)
library(huxtable)
library(nnet)

# Data:
df <- read_dta("Data_for_Replication.dta")

# Set plot theme:
steph_theme <- list(
  hrbrthemes::theme_ipsum_rc(base_size = 12,
                             grid= "",
                             plot_title_size = 20,axis_text_size = 18),
  ggsci::scale_color_npg(),
  ggsci::scale_fill_npg(),
  theme(legend.position = "top",
        plot.caption = element_text(hjust = 0, face= "italic",size = 12.5), #Default is hjust=1
        plot.title.position = "plot", #NEW parameter. Apply for subtitle too.
        plot.caption.position =  "plot"))


source("02. Cleaning-and-recodes.R")

# Materials and Methods / S.I. Balance Checks  -----------------------------------------------------------------------------

# Compared to Control:
regression <- multinom(messages ~ gender + educ + age + pid3 + religpew, data = df, na.rm = T) # Run regression
summary <- summary(regression)$coefficients/summary(regression)$standard.errors #store summary
p <- (1 - pnorm(abs(summary), 0, 1)) * 2 #two tail z test
p

# Compared to Moral Frame:
regression <- multinom(messages2 ~ gender + educ + age + pid3 + religpew, data = df, na.rm = T) # Run regression
summary <- summary(regression)$coefficients/summary(regression)$standard.errors #store summary
p <- (1 - pnorm(abs(summary), 0, 1)) * 2 #two tail z test
p

# Materials and Methods / S.I. Demographics  -----------------------------------------------------------------------------
wpct(df$gender_recode, weight = df$weight_bornagain)
wpct(df$pid_recode, weight = df$weight_bornagain)
wpct(df$educ_recode, weight = df$weight_bornagain)
wpct(df$age_recode, weight = df$weight_bornagain)

# Figure 1 -----------------------------------------------------------------------------------------------------

mod1 <- lm(dx_scale_recode ~ messages * pid7_recode, data = df, weights = weight_bornagain,subset = messages != "Moral")

plot1 <- interplot(m = mod1, var1 = "messages", var2 = "pid7_recode",facet_labs = c("Masks as Helpful"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="Figure 1. Patriotism message effects relative to the control group, by partisanship", 
       x="", y = "Average Effect",subtitle = "", fill = "Gender:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=45,size=7)) +
  scale_x_continuous(breaks = c(1,4,7),
                     label = c("Strong Democrat","Independent","Strong Republican"))

mod2 <- lm(D5_recode2 ~ messages * pid7_recode, data = df, weights = weight_bornagain,subset = messages != "Moral")


plot2 <- interplot(m = mod2, var1 = "messages", var2 = "pid7_recode",facet_labs = c("Masks as Important"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="", y = "Average Effect",subtitle = "", fill = "Gender:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=45,size=7)) +
  scale_x_continuous(breaks = c(1,4,7),
                     label = c("Strong Democrat","Independent","Strong Republican"))


mod3 <- glm(D7_Combine ~ messages * pid7_recode, data = df, weights = weight_bornagain, family = "binomial",subset = messages != "Moral")

# For a smaller output:
# plot3 <- interplot(m = mod3, var1 = "messages", var2 = "pid7_recode",facet_labs = c("Masks Mandate"),
#                    point = T,ci = .90) +
#   steph_theme +
#   theme(axis.text.x  = element_text(angle=90)) +
#   geom_hline(yintercept = 0, linetype = "dashed") + 
#   labs(title="", 
#        x="", y = "Average Effect", caption = "\nMarginal effects of the patriotism treatment relative to the control group by the respondent's partisan identification. 90% confidence intervals (two-tailed) are illustrated. \nThe marginal effects are from a set of analyses in which the treatment conditions were interacted with partisan identification. OLS was used for Masks as Helpful, Masks as Important, and Trump Mask, while Logit was used for Mask Mandate.", fill = "Gender:") + 
#   theme(legend.position = "none",axis.text.x = element_text(angle=45,size=7)) +
#   scale_x_continuous(breaks = c(1,4,7),
#                      label = c("Strong Democrat","Independent","Strong Republican"))

# Larger output:
plot3 <- interplot(m = mod3, var1 = "messages", var2 = "pid7_recode",facet_labs = c("Masks Mandate"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="", y = "Average Effect", caption = "Marginal effects of the patriotism treatment relative to the control group by the respondent's partisan identification. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith partisan identification. OLS was used for Masks as Helpful, Masks as Important, and Trump Mask, while Logit was used for Mask Mandate.", fill = "Gender:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=45,size=7)) +
  scale_x_continuous(breaks = c(1,4,7),
                     label = c("Strong Democrat","Independent","Strong Republican"))


mod4 <- lm(D6_recode2 ~ messages * pid7_recode, data = df, weights = weight_bornagain,subset = messages != "Moral")


plot4 <- interplot(m = mod4, var1 = "messages", var2 = "pid7_recode",facet_labs = c("Trump Mask"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="", y = "Average Effect" , fill = "Gender:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=45,size=7)) +
  scale_x_continuous(breaks = c(1,4,7),
                     label = c("Strong Democrat","Independent","Strong Republican"))


ggarrange(plot1,plot2,plot3,plot4)


# Table 1 -----------------------------------------------------------------------------------------------------

df %>%
  mutate(messages = fct_relevel(df$Exp3_recode, "Control")) -> df

mod1 <- lm(dx_scale_recode ~ messages, data = df, weights = weight_bornagain)
mod2 <- lm(D4_recode2 ~ messages, data = df, weights = weight_bornagain)
mod3 <- lm(D5_recode2 ~ messages, data = df, weights = weight_bornagain)
mod4 <- lm(D6_recode2 ~ messages, data = df, weights = weight_bornagain)
mod5 <- glm(D7_Combine ~ messages, data = df, weights = weight_bornagain, family = "binomial")


huxreg("(Helpful)" = mod1, "(Wear)"  = mod2, "(Important)" = mod3, "(Trump)" = mod4, "(Logit-Mandate)" = mod5, 
       stars = c(`*` = 0.1, `**` = 0.05),error_format = "({std.error})", 
       note = " ** p < 0.1;  * p < 0.05")

# Footnote 7 Table -----------------------------------------------------------------------------------------------------

mod1 <- lm(dx_scale_recode ~ messages * attend_recode, data = df, weights = weight_bornagain)
mod2 <- lm(D5_recode2 ~ messages * attend_recode, data = df, weights = weight_bornagain)
mod3 <- glm(D7_Combine ~ messages * attend_recode, data = df, weights = weight_bornagain, family = "binomial")
mod4 <- lm(D6_recode2 ~ messages * attend_recode, data = df, weights = weight_bornagain)
mod5 <- lm(D4_recode2 ~ messages * attend_recode, data = df, weights = weight_bornagain)

huxtable::huxreg("(Helpful)" = mod1, "(Important)" = mod2, "(Trump)" = mod3, "(Logit-Mandate)" = mod4, "(Wear)" = mod5,
                        stars = c(`*` = 0.1, `**` = 0.05),error_format = "({std.error})", 
                        note = "* < 0.05; ** < 0.01 \n Note: All p-values are two-tailed. Standard errors in parentheses.")


# Footnote 7 Figure -----------------------------------------------------------------------------------------------------

mod1 <- lm(dx_scale_recode ~ messages * attend_recode, data = df, weights = weight_bornagain,subset = messages != "Patriotism")
mod2 <- lm(D5_recode2 ~ messages * attend_recode, data = df, weights = weight_bornagain,subset = messages != "Patriotism")
mod3 <- lm(dx_scale_recode ~ messages * attend_recode, data = df, weights = weight_bornagain,subset = messages != "Moral")

plot1 <- interplot(m = mod1, var1 = "messages", var2 = "attend_recode",facet_labs = c("Religious Message on Masks as Helpful"),
                   point = T,ci = .95) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title = "Average Effect of the Treatment by Attendance", x="", y = "",subtitle = "", fill = "Gender:", "YouGov Survey 2020. 95% Confidence Intervals") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=0,size=7))  +
  scale_x_continuous(breaks = c(1,3,5),
                     label = c("Never","Once a month","More than once a week"))


plot2 <- interplot(m = mod2, var1 = "messages", var2 = "attend_recode", facet_labs = c("Religious Message on Masks as Important"),
                   point = T,ci = .95) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title = " ", x="", y = " " , fill = "Gender:") + 
  theme(legend.position = "none", axis.text.x = element_text(angle=0,size=7)) +
  scale_x_continuous(breaks = c(1,3,5),
                     label = c("Never","Once a month","More than once a week"))


plot3 <- interplot(m = mod1, var1 = "messages", var2 = "attend_recode",facet_labs = c("Patriotism Message on Masks as Helpful"),
                   point = T,ci = .95) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title = " ", x="", y = " ", fill = "Gender:", "") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=0,size=7))  +
  scale_x_continuous(breaks = c(1,3,5),
                     label = c("Never","Once a month","More than once a week"))

ggpubr::ggarrange(plot1,plot2,plot3,nrow = 3, ncol = 1)






